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Abstract 

Binding and unbinding of ligands to specific sites 
of a macromolecule are one of the most elementary 
molecular interactions inside the cell that embody 
the computational processes of biological regula- 
tions. The interaction between transcription fac- 
tors and the operators of genes and that between 
ligands and binding sites of allosteric enzymes are 
typical examples of such molecular interactions. In 
order to obtain the general mathematical frame- 
work of biological regulations, we formulate these 
interactions as finite Markov processes and estab- 
lish a computational theory of regulatory activities 
of macromolecules based mainly on graphical anal- 
ysis of their state transition diagrams. The contri- 
bution is summarized as follows: 

(1) Stochastic interpretation of Michaelis-Menten 
equation is given. 

(2) Notion of probability flow is introduced in rela- 
tion to detailed balance. 

(3) A stochastic analogy of Wegscheider condition 
is given in relation to loops in the state transition 
diagram. 

(4) A simple graphical method of computing the 
regulatory activity in terms of ligands' concentra- 
tions is obtained for Wegscheider cases. 

(5) A general comutational scheme of the station- 
ary probability distribution is obtained in terms of 
probability flows. 

keywords: Biological regulation, molecular interac- 
tion, Markov process, stationary probability distri- 
bution, probability flow 



1 Introduction 

Control is a ubiquitous built-in mechanism sup- 
porting a variety of functions of living organisms. 
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The advent of cybernetics in the late 1940's estab- 
lished a close link between control in living organ- 
isms and that of man-made artifact. Notions like 
feedback, feedforward, robustness, stability, and 
noise, which are fundamental in control engineer- 
ing turned out to be relevant to the control mecha- 
nisms of living organisms. Rapid progress in molec- 
ular biology in the last fifty years opened up a new 
world of intracellular control of biochemical pro- 
cesses including those of metabolism and gene ex- 
pression. Studying control mechanisms became an 
important issue of biology. Monod correctly called 
this field microscopic cybernetics [31j . while others 
called it regulatory biology [50] . The continuing en- 
deavor to clarify the material basis of regulations 
reveals vast complexity of regulatory networks as- 
sociated with huge variety of material links inside 
a cell, which in turn motivated considerable inter- 
est in the quantitative analysis of regulation mech- 
anisms through mathematical modeling and sim- 
ulation. The modeling/simulation paradigm has 
been used to reveal the underlying structural prop- 
erties of complex regulatory networks. Existence of 
cellular switches [JJ [8] [M] [23] [30] [34] [37] , evaluation 
of molecular fluctuations [TH] [H] [H] [13] > analysis of 
biological oscilations [5] [7] [17] , consideration of net- 
work robustness [35] [5T] , theoretical demonstration 
of cellular behaviors [TTJ [T2] [35] , are examples of its 
contributions. Excellent reviews of these achieve- 
ments are found in [UJ [M] [H] ■ Recently, we pro- 
posed the notion of compound control [42] which 
captures a salient characteristic feature of biologi- 
cal control. Compound control is a biological way 
of realizing proper input/output relationship of reg- 
ulatory dynamics to cope with diverse, complicated 
and unpredictable environmental changes. 

It is paradoxical that the theoretical and/or phys- 
ical basis of the modeling/simulation paradigm has 
not been well exploited, in spite of the remark- 
able success it has attained. For instance, in many 
papers on intracellular regulations, the traditional 
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Michaelis-Menten-Hill equation in enzymology is 
extensively used for describing molecular interac- 
tions, e.g., the action of repressors in transcrip- 
tion regulation, or the effect of feedback [5] [2] [57] 
L 30j|34 48J, without sound justification. 

The transcription rate is regulated through bind- 
ing of transcription factors in some domains of 
DNA, the so-called cis-regulatory elements [6], that 
are able to affect transcription initiation. The affin- 
ity of cis-regulatory elements to the transcription 
factors determines the overall transcription rate of 
the gene. They are correlated with each other and 
collectively exhibit synergistic effects through which 
a complex computation is performed to achieve 
adequate transcription control [26] . These fac- 
tors are represented in the context of a thermo- 
dynamical equilibrium energy distribution suggest- 
ing to plausible models of transcription regulation 
[3] [4] [35] [47] . This paper aims to establish a unified 
framework with solid mathematical grounding and 
physical plausibility that is able to capture various 
bio-physical attributes of intracellular regulations. 

A significant portion of genes encode enzymes 
which are involved in metabolic control, another im- 
portant form of intracellular regulation. Metabolic 
control is a relatively old subject of bio-chemistry, 
where allosteric enzymes are the major agents of 
control. Allosteric enzyme is a macro-molecule 
that has several binding sites which can bind sub- 
strate and ligands. The binding pattern of ligands 
changes the conformations of the enzyme, through 
which the activity of the enzyme for the target re- 
action is controlled. Quantification of the activity 
changes of allosteric enzyme has been a central issue 
in biochemistry since the classical Monod-Wyman- 
Changeux model was first proposed [32] [49]. Since 
then, various theoretical and experimental methods 
have been used to predict the conformation changes 
due to substrate binding and modification by lig- 
ands [3] [9] [T3] [T5] [IE] [24] [29] [36] [40] [46] [52] . 

The regulatory mechanism of allosteric enzymes 
is in some sense similar to that of operons, in spite 
of the obvious functional and structural differences 
between them. Both of them are controlled by in- 
teractions between binding sites (catalytic sites for 
metabolic cases and cis-regulatory elements for ge- 
netic cases) and binding factors (ligands and sub- 
strate in the metabolic case and transcription fac- 
tors in the genetic case). In the metabolic case, in- 
teractions take place between proteins, while in the 
genetic case they take place between protein and 
DNA. It is interesting that although the allosteric 
protein is a player in both regulations, it plays to- 
tally different roles. In the metabolic control, it is 



controlled by ligands, while in the genetic control it 
acts as transcription factors to control genes. 

The protein/protein interactions of metabolic 
regulation and protein/DNA interactions in tran- 
scription control are weak and essentially reversible, 
and they share many properties. In fact, equi- 
librium statistical mechanics is used to quantify 
both metabolic regulation [15] and genetic regula- 
tion [4] [35] [47]. If we view the binding/unbinding 
processes taking place between binding sites and 
binding factors in both cases as a common stochas- 
tic process, it suggests a common theoretical frame- 
work to deal with metabolic and genetic regulations 
in a unified way. 

Along this line, we propose a common theo- 
retical framework for describing biological regula- 
tors of metabolic and genetic regulations to faci- 
litate quantitative descriptions of intracellular bio- 
chemical processes. The core of our approach is 
to describe protein/protein and protein/DNA in- 
teractions with a common Markov process focusing 
on the behavior of a single molecule rather than a 
population of molecules. This viewpoint enables us 
to describe phenomena with a finite state Markov 
process, which greatly simplifies the argument by 
avoiding the difficulties associated with infinite di- 
mensionality [33]. All the quantification features 
associated with various complicated physical and 
biological phenomena can then be condensed into 
transition probabilities. We introduce a new notion 
of probability flow that offers a new insight into the 
stationary distribution of the Markov process. The 
stationary equation is regarded as representing the 
conservation of probability flows at each node. It 
is closely related to the Wegscheider condition de- 
rived a century ago. Based on the notion of prob- 
ability flow, we derive a representation of the sta- 
tionary probability distribution, which leads to a 
unified quantitative form for the general biological 
regulator. Our method of computing the station- 
ary probability distribution dramatically simplifies 
the classical King and Altman method [53], which 
is used for computing the stationary distribution of 
chemical processes. It is expected to yield a new 
computational tool for operon regulation that can 
deal with the increasing complexity of operons |38] . 

Metabolic and genetic controls have been inves- 
tigated relatively independently because of their 
apparent large differences in modality. Therefore, 
there have been few studies dealing with systems 
that combine metabolic and genetic regulations 
[52] , The unified framework described in this pa- 
per will enhance our understanding of intracellular 
biochemical processes. 
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In the next section, we formulate a biological 
regulator in an abstract way as a Markov process. 
The abstract biological regulator is defined through 
the stationary solution of the Master Equation. In 
Section [3l we introduce the notion of probability 
flow and explain its relevance to the stationary dis- 
tribution. Sections 0] and [5] deal with the single 
loop case and introduce the Wegscheider condition, 
which guarantees that the probability flow vanishes. 
Section [6] generalizes the results of the preceding 
section to multi-loop cases. Section [7] discusses the 
meaning of transition probabilities to deduce bi- 
ologically meaningful representations of biological 
regulations. 

2 Characterization of biologi- 
cal regulations as a finite- 
state Markov process 

The main stage of genetic control is transcrip- 
tion regulation whose fundamental mechanism is to 
change the transcription rate of the operon via bind- 
ing/unbinding of transcription factors to and from 
their cis-regulatory elements. The main agent of 
metabolic control is the allosteric enzyme that de- 
termine the rate of the target chemical reaction via 
binding/unbinding of substrate and ligands. The 
fundamental common feature of both regulations 
lies in the molecular interaction between sites of the 
macro-molecule and molecular binding factors that 
work through binding/unbinding processes (see Ta- 
ble [J). 

These processes are obviously random subject to 
certain thermodynamical constraints. Our idea is 
to establish a mathematical framework to describe 
such regulatory mechanisms in a unified way. We 
now formulate the biological regulation in a some- 
what abstract way. We assume that the macro- 
molecule, which is the main agent of the regula- 
tion, has n binding sites (BS), each of which can 
bind some of m ligands with different affinities that 
depend on the binding patterns of sites. Instead of 
the ligand, we use the word binding factors (BF) 
to emphasize the bilateral symmetry of sites to be 
bound and factors to bind. Each BS can bind only 
one BF at each time. The binding sites are de- 
noted by 61,627- •• ,b n , while the binding factors 
by U\,U2,--' 7 U m . The state of the regulator is 
denoted by n-tuples of m alphabet Ui, U^, ■ ■ ■ , U m 
denoting the BFs being bound and <j> which repre- 
sents the empty (unoccupied) site. If the regulator 
has three BSs, for example, S = denotes 
the state with b\ and 63 being bound by U\ and U2, 



respectively, and with 62 empty. Let N be the total 
number of non-empty states, which is obviously fi- 
nite. Since we have the empty state So = (4>4> ••■</>), 
the total number of states is N + 1. 

The state of the regulator changes as its bind- 
ing pattern changes, i.e., as BFs bind to and dis- 
sociate from BSs. These are clearly stochastic phe- 
nomena and are legitimately described as stochas- 
tic process. Let the transition probability from the 
state Sj to Si during an infinitesimal time duration 
At be denoted by g^At. All the quantitative fea- 
tures of the complex biochemical process associated 
with protein/protein interactions and protein/DNA 
interactions including origomerization [36 49J , con- 
formational change [13] [35) and target localization 
[9] can be adequately represented in terms of qij. 
For the sake of clarity, we assume that binding or 
unbinding of only one BF to and from a BS occurs 
during the infinitesimal time duration At. Thus, 
the transition from {Ui^U^) to (UiipUs) is regarded 
as consequtive transitions from (Z/i^f/a) to (Ui4>4>) 
and from (Ui<fxf>) to (E/i^L/3). The transition prob- 
ability for unit time represents the rate of reac- 
tion from Sj to Si, and sometimes it is regarded as 
identical to the rate constant of the corresponding 
reaction. We define Aj to be the set of states that 
are accessible from Si during the infinitesimal time 
duration At. We call it the adjacent set of S{. 

Let pi(t) be the probability that the regulator 
is in a state Si at time t. Then, the stochastic 
time-evolution of the regulator is described by the 
chemical master equation (CME), 

-J: = H QijPj ~ H 9jiPi- i = 0, l, • • • ,N. (1) 

The notation X^jeA denotes the sum of all j such 
that Sj £ A,. The first term of the right-hand side 
of {1} represents the net probability coming from 
adjacent states to Si while the second term the net 
probability leaving Si to adjacent states. The CME 
([I]) is written in the matrix form as 

dJ dr = W) (2) 

where pit) denotes the (N + l)-dimensional vector 
whose (i + l)-th component is Pi(t). The matrix Q 
in ([2]) is called a transition matrix. 

A finite Markov process can be described by the 
state-transition diagram. In Figure [T] shows exam- 
ples of state-transition diagrams. The most salient 
feature of the Markov process for our purpose is 
the reflection property that q^ ^ always implies 
qji 7^ because of the reversibility of the molecu- 
lar interactions we are interested in this paper. In 
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Table 1: Comparison of genetic and metabolic regulations. 





Agent 


Site 


Factor 


Control 
Parameters 


Genetic 


Operons 


Cis-regulatory 
elements 


Transcription 
factors 


Transcription 
rate 


Metabolic 


Allosteric 
enzymes 


Catalytic sites 


Ligands and 
substrates 


Reaction rate 




Figure 1: State transition diagram of (a) Example 
1, (b) Example 2,(c) Example 3. 



terms of the adjacent set A,;, Sj £ A; always im- 
plies Si £ Aj. 

It is known [JS] that the solution of © converges 
to a unique equilibrium stationary solution p that 
satisfies 

Qp = o, (3) 



or equivalently, 



Y QjiPi 



(4) 



Equaiton l[3" |) or is called the stationary state 
equation (SSE). Usually, the convergence is faster 
than the chemical reactions controlled by the regu- 
lator. Therefore, we can assume that the regulator 
is always in a stationary distribution satisfying (|3j). 
together with the normalization constraint, 



N 



5> = 1 - 



(5) 



Let 7j be the activity of the regulator at state 
Si. The overall activity 7 of the regulator is then 
defined as the average activity with respect to the 



stationary distribution p, i.e. 



N 
i=0 



(6) 



The examples that follow show that the above def- 
inition of the regulation activity adequately de- 
scribes both genetic and metabolic regulations. 

Example 1. (Stochastic version of Michaelis- 
Menten equation) 

Consider a molecule of an enzyme E that cataly- 
ses the reaction of producing a product P from a 
substrate S, i.e., 



ki 

E + S ES 

k-i 



E + P. 



(7) 



An enzyme molecule can be regarded as a biolog- 
ical regulator with one BS and one BF (the sub- 
strate). The enzyme has the two state Sq and Si, 
the empty state and the state occupied by a sub- 
strate molecule, respectively. The transition dia- 
gram is shown in Fig. [TJa). The transition matrix 
Q in this case is written as 



Q 



-qw 
qio 



-<7oi 



and the stationary distribution is given by 
1 



Po 



1 



£10 

<3oi 



Pi 



'/HI 



1 



301 



(8) 



(9) 



Since (?oi represents the probability of unbinding the 
substrate S from the BS during unit time, we can 
identy it with k—\. In addition, the probability of 
binding S to the BS is proportional to the concen- 
tration [S] of the substrate with ki being its rate 
coefficient. Thus, 



qoi 



k-i, qi = k 1 [S], 



(10) 



which gives 



Po 



K x + [S] 



Pi 



[S] 



K x + [S] ■ 



ki 



(11) 
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Since pi denotes the probability of the state of the 
enzyme with S being bound (ES), the total number 
of ES molecules is given by p\ [E] T , where [E) T de- 
notes the concentration of the total enzyme, which 
is assumed to be constant. Therefore, since the 
maximum reaction rate is given by k 2 [E] T , the ra- 
tio of the reaction rate v to its maximum v max is 
given by 



k 2 p\ [gjj 
k 2 [E] T 



= Px = 



[S] 



K x + [S] ' 



(12) 



which is the celebrated Michaelis-Menten equation. 
Here, the rapid equilibrium assumption used in de- 
riving the Michaelis-Menten equation [3U] has been 
replaced with the notion of rapid convergence of 
the probability distribution to a stationary one. 
Note that we consider the behavior of a single en- 
zyme molecule rather than the collective behav- 
ior of the enzyme and enzyme-substrate complex. 
In this way, a stochastic formulation of biologi- 
cal regulation gives an alternative interpretation of 
the Michaelis-Menten equation. It should be noted 
that we have not assumed the enzyme concentration 
is constant, which is required to derive Michaelis- 
Menten equation in the traditional way. 

It is not difficult t o see that the standard devia- 
tion a v = \J v 1 — v 2 is given by 



K x + [S] ' 



which gives a rough estimate of the precision of the 
approximation (fT2"|) . This is a bonus of the stochas- 
tic version of the Michaelis-Menten equation that 
we have derived. 

Example 2. (Allosteric enzyme with an activator 

m) 

An allosteric enzyme e with one binding site for an 
activator A in addition to the one for a substrate S 
is described as 

E + S i± ES h E + P 
fe-i 



+ 
A 

k- 2 ][k 2 

EA + S 



+ 
A 

k-sf-Ua 



(13) 



ESA k ^ EA + P 



The enzyme has now the two binding sites, one for 
the substrate and the other for the activator A. The 
enzyme has four states So = (4>4>), Si = (S<fi), S2 — 
(SA), S3 = (4>A). The state transition diagram is 



shown in Fig. [ljb). We can associate the transition 
probabilities with the rate constants in (fT3|) as 



qxo = fa [S] , q x = k-x, 921 

932 = fe-4, 923 = ^ [S] , 903 



^3 [-4] , 912 = k- 3 , 
k-2, 930 = k 2 [A] , 
(14) 



where [S] and [A] denote the concentrations of the 
substrate and the activator, respectively. The tran- 
sition matrix is given by 



Q 



D 
9io 


930 



9oi 
Dx 

921 







912 

D 2 

932 



903 



923 

D 3 



(15) 



where D = -(qxo + 930), D 1 = -(q x + 921), D 2 = 
-(912 + 932), D 3 = -(903 + 923)- The rate of the 
target chemical reactions producing P is given by 

7 = kpPi + k pA p 2 

The actual computation of the stationary probabil- 
ity distribution for this example is done in Example 
6, where you see the result is very complex. 

Example 3. (Allosteric enzyme with both activa- 
tor and inhibitor) 
The reaction scheme is written as 



EI + P 



EI - 


f S 


^ ESI ^ 

fe-7 




6 


fetlk-s 


I 




I 


+ 




+ 


E + 


S 


fei u 

£ ES ^ 
fc-1 


+ 




+ 


A 




A 


k- 2 Uk 2 






EA 


f S 


4 ESA ^ 

k—4 



E 



P 



(16) 



EA + P 



where I denotes the inhibitor. The enzyme has six 
states, namely, So — (<t>4>), S\ — [4>A), S 2 — (SA), 
S 3 = (S<f>), Si = (SI), S 5 = (4>I). The state transi- 
tion diagram is shown in Fig. []Jc). The right part 
is identical to the diagram of Fig. [lja). The state 
transition matrix is given by 



Q 



'Do 


9oi 





903 





905 


9io 


Dx 


912 














921 


D, 


923 








930 





932 


D 3 


934 














943 


Di 


945 


950 











954 


D 5 



(17) 
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where D = -{qio + q 30 + qso), D x = -{qoi + q 2 i), 
D 2 = ~(qi2 + 932), D 3 = -(q 3 + q 23 + 943), D 4 = 
-(934 + 954), D 5 = -(g 5 + 945)- The rate of the 
target chemical reactions is given by 

7 = kp A P2 + kpp 3 + k pI p 4 

The computation of the stationary probability dis- 
tribution is done in Example 6. 

Example 4. (Lac operon) 

The lac operon, which has been a subject of ge- 
netic molecular biology for more than fifty years 
([5] [38] [48] [52]), can be regarded as a biological reg- 
ulator. A number of different models have been 
proposed for it including a very complicated model 
with seven cis-regulatory elements |38j . The stan- 
dard model has three binding sites [55]. One is the 
promoter with RNAP as its unique binding factor. 

There are two binding factors associated with the 
lac operon, cAMP-CRP and LacI, which have their 
own binding sites. They are independent of each 
other. cAMP-CRP is an activator of the operon, 
while LacI is a repressor. Binding of LacI at its 
site prevents the other factors from binding to their 
sites. 

Denote the RNAP, cAMP-CRP and LacI by U x , 
U 2 and U 3 , respectively. Then, the lac operon has 
the following five states; 

So = (#0), Si = (UkM), S 2 = (UiUifi, 

s 3 = (0^20), s 4 = (<m>u 3 ). 

The transition diagram is shown in Fig. O The 
transcription rate is given by 

7 = 7iPi + 72P2, 

where 71 and 72 are the transcription rates corre- 
sponding to Si and S 2 . 




Figure 2: State transition diagram of lac operon. 

Recently, Santillan et al. proposed a lac operon 
model with five binding sites [38j . The number of 
binding patterns of its cis-regulatory units is at 
least 50. Figure [3] shows the transition diagram 
which is very complicated. We need an efficient 
computational theory to handle complex operons 
like this. 




Figure 3: Transition diagram of lac operon with 50 
states. 

3 Conservation of probability 
flows in stationary distribu- 
tion and loop-free transition 
diagram 

In this section, we introduce the notion of probabil- 
ity flow and discuss its meaning in solving the SSE 
QJ. The SSE (0]) can be written as 

^2 (QvPj - QjiPi) = °> i = 0, 1, • • • , JV. (18) 

The term qijPj — qjiPi represents the net probability 
of the state transition from Sj £ Aj to Si. Hence, 
it is reasonable to call it the probability flow from 
Sj to Si, which is denoted by 

Pij = qijPj - qjiPi- (19) 

If > 0, the probability flows from Sj to Sj and 
if < 0, it flows oppositely. Clearly, it is skew- 
symmetric, i.e., 

Pij + pji - 0. (20) 

The SSE ([P8]) is represented in terms of probability 
flow as 

E^ = ' Vi ' ( 21 ) 

which implies that the probability flows are con- 
served at each state node. In other words, the net 
probability flow incoming to Sj (the sum of positive 
Pij) is equal to the outgoing flow from Sj (the sum 
of negative pij), as is illustrated in Figured) Here, 



G 




Figure 4: Conservation of probability flows. 



it is important to notice that the probability flow is 
directed. 

In chemical kinetics, one usually assumes that the 
equilibrium state satisfies the relations 



qijPj = <i,,i>,. Vj € At, Vi, 



(22) 



if we interpret as the kinetic rate coefficient of 



the reaction Si 



Si. The relations 



are called 



detailed balance in chemical kinetics |I9j . Due to 
(IT5|) . the detailed balance holds if and only if = 
0, Vi, j G A,; i.e., all the probability flows vanish. 
The relation (1221) is written as 



The direction of an edge can be freely assigned in 
TR diagrams, but it must be consistent with the 
direction of the TR, in the sense that an edge with 
rij as its TR must be directed from Sj to $£. Now, 
we shall prove that the probability flows vanish at 
all edges unless the transition diagram has a loop. 
Here, a loop is defined in the context of the TR 
diagram. In other words, if the transition diagram 
does not have a loop (loop- free), all the probability 
flows vanish. To see this remarkable fact, assume 
that there exists an edge e with a non-zero proba- 
bility flow and a state Sj is connected to this edge. 
Then, due to the conservation of probability flow 
at Sj, there exists another edge connecting Sj to 
another state Sj with a non-zero probability flow. 
The same reasoning for Si leads one to conclude 
that Si must be connected to a new state Sk by an 
edge with a non-zero probability flow (Fig. [6|) . Re- 
peating the procedure creates a sequence of states 
Sj — ► Si — ► Sk — > • • • , which are connected by edges 
with non-zero probability flows. Since the number 
of states is finite, the sequence must visit a state 
which has already appeared in the sequence. Hence, 
the graph must contain a loop. 



Pi 



r ijPj 



(23) 



where m is called the transition ratio from Sj to 
Si and is defined as 



— -is, 

e 



Si 



Oil 



(24) 



Figure 6: Sequence of edges with non-zero proba- 
bility flow. 

For a while, we concentrate on the loop-free case 
where all the probability flows vanish. Take two 



The transition ratio (TR) is directed as shown in 
Fig. E^a) and corresponds to the equilibrium coeffi- 
cient of the chemical reaction S 3 ±5 S z . It is some- arb itrary states Si and s]. If there are two different 
times convenient to describe the transition diagram pathg connecting Sl and Sj, it means that there is 
in terms of the transition ratios (TR), instead of 
transition probabilities, as is shown in Fig. EJb). 
We call such a diagram a TR diagram to distin- 
guish it from the usual transition diagram. Note 
that 

rv=r~ 1 . (25) 



a loop composed of a path from Sj to Si and one 
from Si to Sj. Therefore, if the transition diagram 
is loop-free, the state Si is connected to Sj through 
a unique path. Let this path be composed of I state: 



Si t — Si 



Sj — Si 1 

the transition probability from Sj to Si as 



We can define 



la) 



Figure 5: TR Diagrams. 
TR diagram of Fig. QTb) 




(a) Description of TR, (b) 



Hi 



(26) 



where the overbar is used to distinguish the transi- 
tion probability computed along a path from those 
between adjacent nodes. In this way, we can com- 
pute the transition probability for any state pair Si 
and Sj which are not necessarily adjacent to each 
other. 

We can extend the detailed balance ((22)) to any 
pair of states by using the extended transition prob- 
abilities (1261). 
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To see this, assume that three states So, Si, and 
S2 are connected by a path, i.e., S\ G Ao and 
52 G Ai. Then, the detailed balance (|22p implies 
QioPo = qoiPi and g 2 iPi = <?12P2- Multiplying each 
side yields qioPo^iPi = <?oiPi<7i2P2- Cancelling p 1 
from both sides gives q2iqwPo = Q01912P2, or equiv- 
alently, <?2oPo = qmPi- This relation suggests that 
the detailed balance holds in terms of the general- 
ized transition probability (|26[) even for state pairs 
which are not adjacent to each other. It is not diffi- 
cult to see that this is indeed the case by repeatedly 
applying it through the unique path connecting the 
two states, i.e., 



QijPj = QjiPu Vi,j. 



(27) 



The above relation includes the detailed balance 
(|22|) as a special case. Therefore, we call the re- 
lation (|27| the generalized detailed balance. Taking 
j = in (f2"7) , we can represent pj as 

Pi = r i0 p , (28) 

where fy is defined as a generalization of (|24p with 
qij being replaced by q~ij given by (|26p . i.e., 



q~ji 



(29) 



Here, the overbar is again used to distinguish it 
from Tij given in (|24|) for state pairs adjacent to 
each other, fij is also called the transition ratio 
(TR) from Sj to Si. Note that the relation (|25|) is 
extended to 

r ij = r ji ■ 

From the normalization constraint and (j2"5)) , the 
statitonary probability distribution is simply given 

by 

Pi = y , i = 0,l,-'-,N (30) 
E J= o r jo 

where foo = 1- 

We sum up the discussion in this section as fol- 
lows: 

If the transition diagram is loop-free, the following 
facts hold: 

(1) the probability flows vanish at every edge. 

(2) The stationary distribution is given by H3U\) . 

(3) The generalized detailed balance holds for 
the stationary distribution. 

Example 5. (Loop- free Transition Diagram) 
Consider a transition diagram of Fig. [7j which is 
loop- free. According to (f28|) . we have p\ — qiopo, 
P2 = r 2 ir w p , P3 = r 3 ir w po, Pa = riz^ipo, P5 = 
7"5oPo with 

1 

P ° 1 + rio(l + r 2 i + r 3i + r 43 r 3 i) + r 50 



Figure 7: A loop-free transition diagram. 



The above reasoning can be applied to a loop- free 
subgraph of the general (not loop-free) transition 
diagram. If a loop-free subgraph is attached to a 
state Si which belongs to a loop, we can write down 
the probability of the state of that subgraph ac- 
cording to (j30|) with p being replaced by pi. As an 
example, consider the transition diagram of Fig. [8l 
The state Sj is a part of a loop which will be dis- 
cussed in the next section. If the state Si is not in 
any loop, then it is a part of a path connecting Si 
and the terminal state Sfe. Thus, 



Pi 



jPji 



(31) 



because the probability flow at any edge on the path 
connecting Sj and Sk vanishes. This can be shown 
by directly applying the argument of this section to 
the subgraph of Fig. [H 




Figure 8: A loop-free subgraph. 



4 Wegscheider condition 

In the preceding section, we showed that the proba- 
bility flows vanish if the transition diagram is loop- 
free. In that case, the stationary distribution is 
simply calculated from ([50)) . There are cases where 
the probability flows vanish even if the transition 
diagram contains loops. 

In order to discuss this issue, we consider the case 
where the whole transition diagram is a loop, as 
shown in Fig.[9ja). The corresponding TR diagram 
is shown in Fig. (9^b). The state transition matrix 
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Figure 9: Loop transition Diagram, (a) Transition 
probability description, (b) transition ratio descrip- 
tion. 

is given by 





9oi 


• 





qoN 


9io 


D x 


9l2 • 











921 


D 2 ■ 













• 


QN,N-1 


D N 



with Di = — + qi-i,i), i = 1,2, ■ • • ,N — 1, 
-Do = -(gio + qm), D N = -(qoN + gjv-i, AO- 
Each state has only two edges, one of which re- 
ceives the probability flow, the other dispatches it. 
Due to the conservation of probability flows (f2Tj) . 
the incoming probability flow and the outgoing one 
are equal. Hence, each edge of the loop has the 
same amount of probability flow, which is denoted 
by p. Taking its direction to be clockwise, that is, 
taking the probability flow from Sj_i to Si to be 
positive, the probability flow is given by 

P = qis-iPi-i ~ q%-\,iPi, i = 1,2, ■•• , N. (33) 
Thus, we have the following recursion: 

Vi = ri,i-iPi-i — , i = l,2,---,N, (34) 

Qi-l,i 

P 

Pa = roNPN ■ 

qN,o 

In order to obtain a general representation of 
stationary probabilities, we introduce the numbers 
£o, £i, ' " ■ j £n defined sequentially by 

£o = o 

Si = ^-1^-1 — , i=l,2,.-.,JV (35) 

Qi-l,i 



Then, due to ((31]) and ((23), Pi = r^-iPi-i + p(£i - 
r M-i £2-1)1 Therefore, we have 

Pi ~ p£,i = ^.i-iOi-i - p£i-i), 

which implies p l - p& = f l0 p . Here, f i0 = 
ri,i-iTi-i,»-2 ■ • • rio- Thus, pi is represented as 

Pi — fiopo + p^, i=l, ■■■ ,N. (36) 

If p = 0, the above relations are identical to 
(|28[) . Thus, equation (f3"6"| represents the stationary 
distribution as a sum of the probability distribu- 
tion for the case with zero probability flow and a 
correction term due to non-zero probability flow. 
Since p — qoNPN — 9jvoPoj the relation (|36|) for 
i = N implies p = q QN {fNoPo + p£n) - qmPo = 
qNo(r ON fNo - l)Po + pqoN^N- Therefore, we get 

pO- - qoN^N) = -(1 - v)qNopo (37) 

where 77 is given by 

V = r 0N f N0 = r 0N r N . N -x ■ ■ ■ r w (38) 

Since £1 = — I/901 < and r.y > 0, £at < 0. Hence, 
1 — 5oiv£zv > 0. Therefore, we conclude that the 
probability flow vanishes, i.e., p — 0, if and only if 
77 = 1. 

The number 77 is the product of all the TRs along 
the loop. If we use the analogy between the TR and 
the equilibrium coefficient in chemical kinetics, the 
condition 77 = 1, i.e., 

rior2i • • • r N = 1, (39) 

corresponds to the condition derived almost a cen- 
tury ago by Wegscheider. The condition (|3"9")) is re- 
ferred to as the Wegscheider condition [22 . We 
call the product i] of (|38|) the Wegscheider prod- 
uct for a loop as in Fig. [9] The Wegscheider prod- 
uct is in some sense directed. The representation 
([38]) defines the Wegscheider product clockwise in 
the context of Fig. [HI It can also be defined in the 
opposite direction 77' = TnoTn-i^n ■ ■ - rox- Due to 
([25jl. 77' = 77 _1 . Since n = 1 implies ?/ = 1, the 
Wegscheider condition can be represented by com- 
puting the product in either clockwise or counter 
clockwise direction. 

The Wegscheider condition is a property of a loop 
that requires the transition ratio along the total 
loop to be unity. This condition is equivalently writ- 
ten in terms of transition probabilities as 

gio?2i • • ■ qoN = qmqN-i.N • ■ ■ 901 ■ (40) 

The left-hand side represents the transition proba- 
bility from a state to itself along a clockwise con- 
tour, while the right-hand side represents that along 
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a counter clockwise contour. The identity (|4"0"1) gives 
a stochastic interpretation of the Wegscheider con- 
dition. 

Taking the normalization condition ([5]) into ac- 
count, we have 



rio + Pii 



z = 0,l, ••• ,N, (41) 



where roo = 1, Co = 0. If the Wegscheider condition 
holds, then p = and the stationary distribution 
becomes identical to (p0| . 

The most common method of quantifying the 
regulatory activities of operons is based on ther- 
mal equilibrium theory 3J |4J , which assigns a Gibbs 
free energy to each binding pattern of the binding 
sites. This corresponds to assigning a stationary 
probability distribution a priori, rather than con- 
structing a Markov process by assigning the tran- 
sition probabilities between states. The transition 
ratio may then be defined as the ratio between the 
state probabilities. In that case, the Wegscheider 
condition obviously holds. To see this, consider a 
loop So -> Si -> 52 -> S - Then, r w = pi/po, 
T21 = Vi/Pl, r 2 = Po/P2- Thus, r w r 2 iro2 = 1. 
This proves that the Wegscheider condition is con- 
sistent with thermal equilibrium theory. The rela- 
tion ([4Tj) . however, suggests more general types of 
equilibrium. 

5 Edge removal and modifier 

In this section, we consider the meaning of the num- 
bers fa in (|35p and try to interpret ([36)) in the 
context of edge removal. From ([35]) , <7i_i,i£i — 

<?M-l£i-l = — 1 arL d — = ~ 1- 

These relations yield - + 

+ qi-i,i)£i = °- Since Co = and £1 = 
— 1/goii we have, from 



oc = 



-1 



(1 - Qon^n), 



(42) 



where £ — [Co £1 • ■ • £zvl • To further inves- 
tigate the meaning of the vector £, we consider a 
transition diagram which is obtained from the orig- 
inal diagram of Fig. [9] by eliminating the edge con- 
necting So and Sn (Fig. HUT a)). The modified tran- 
sition diagram has no loop, and hence, the station- 
ary distribution is given by pi = fjoPOj following the 
discussion in the preceding section. The transition 



matrix Q' corresponding to the modified diagram 
of Fig. [TOTa) is obtained by taking qoN = <Zjvo = 
in ([31]). From (gSJ), we have 



-i 

o 



1 



(43) 



We call the vector £ satisfying (|43|) with £o = 
the modifier of the loop corresponding to the edge 
Sn — > Sq. It is computed through the recursion 
formula (f35|) . We shall generalize it in the sequel. 




Figure 10: Reduced transition diagrams, (a) Edge 
Sn — > Sq is eliminated, (b) edge S^ — > S v +i is 
eliminated. 

If we let p' be the vector whose (j + l)th element 
Pi is given by (|28p . which represents the stationary 
probability distribution corresponding to the mod- 
ified loop-free diagram of Fig. [TOTa). equation (|36| 
can be written as 



P = P' + Pt 



(44) 



This is the general form of the stationary probabil- 
ity distribution. It holds for a loop transition dia- 
gram as well as for general transition diagram with 
multiple loops, as is discussed in the next section. 
Note that the stationary distribution is represented 
by a stationary distribution p' corresponding to the 
modified loop-free diagram with a correcting term 
p£ due to non-zero probability flow. 

Equation (|44[) is obtained by cutting the edge 
Sn — > So as in Fig. [TOT a). Here, we can show 
that equation (|44|) is also derived by eliminating 
any edge S u — > S v +i of the loop (Fig. flOTb)). Since 
fio is defined as the product of each TR along the 
path connecting Si to Sq, we have 



no 



fio?*2i • • ■ r,-,i_i if i < v. 

rmrN-i.N ■ ■ ■ n,i-i if i>v + l. 



(45) 
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Now, the modifier is defined as 

[Forward Recursion] 
£o = 

& = ri,i_i&_i , i 

Qi-l,i 

[Backward Recursion] 
6v+i = 0. qN+l.N = Oqn 

£i = n,i+i£i+i H , i 

Qi+l,i 



1,2,- 



(46) 



iV,iV-l,--- + 

(47) 



The recursion formula (]46[) computes the modifier 
elements of the path connecting So to SV, while the 
recursion formula (|47|) computes those of the path 
connecting Sq to S^+i. The paths are oppositely di- 
rected, but both recursions are essentially the same 
except the signs of the added term. We call (|4"6"1) for- 
ward recursion, while (|47p backward recursion. It is 

straightforward to see that £ = [£ £1 • ■ • £at] 
satisfies 









Qv+l.v 











Iv.v+l 















" " 













1 






+ 


-1 









_ _ 



Since q u +i.v = Q_v,v+i = in the modified transition 
diagram, we have 




1 

-1 




•v + 1 
■u + 2 



(48) 



It is not difficult to see that the stationary distri- 
bution is given by (|36[) with fjo and £j being repre- 
sented by (|4"5")) and (j4*B"|) (|4"7|l . respectively, by apply- 
ing the arguments of the preceding section to the 
forward and backward recursions separately. 

The modifier is extended to any edge S v —* S v +\ 
and can be computed through (|46p(|T7|) . The sta- 
tionary probability distribution is given in this case 
by (|44p . where p' corresponds to the stationary 
probability distribution of the reduced loop-free di- 
agram. Since p — holds under Wegscheider con- 
dition, we obtain the following result: 



The stationary distribution of a loop is unchanged 
if an edge is removed, provided that the Wegschei- 
der condition i39\) or J^0| ) holds. 
In other words, the computation of the stationary 
distribution for a loop transition diagram is reduced 
to the loop-free case where a very simple form (|3T)1) 
is already available by eliminating an edge, provided 
that the Wegsheider condition holds. We shall ex- 
tend this remarkable property to the general tran- 
sition diagram with multiple loops in the next sec- 
tion. 

Example 6. (Stationary Distribution of Example 
2) 

This is the case of N = 3 in the above algorithm. 
Remove the edge e : S2 — ► Sa(v = 2). The numbers 
£i,i = 0,1,2,3, are given respectively by forward 
recursion 

£o=0, & = - — , 

ttx 1_ 

9oi 912 
and by backward recursion 



6 = - 



921 + 9oi 
9oi9i2 



(49) 



£3 = 



1 

903 



The stationary probability distribution is calculated 
to be 



Pi = r w p 



P2 = r 2 oPo 



1 

— / 

9oi 

1 



r2ir w po 



9oi + 921 
901912 



(50) 



P3 = r 30 Po H P- 

903 



The Wegscheider product 7/ is given by 

903932921910 

f] = r a3 r 32 r2ir w = . 

9309239i29oi 

The probability flow p = 932P2 — 923P3 is given by 

903932921910 — 930901912923 



P = 



901912(903 + 932) + 903923(901 + 921) 



Po 
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Eliminating loops from 
the transition diagram and 
derivation of the flow equa- 
tion 



Now, we consider the general case where multiple 
loops exist, and derive a method of computing prob- 
ability flows. If some edges are removed, the tran- 
sition diagram becomes loop-free. Although the se- 
lection of such loops is not unique, the minimum 
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number of such edges required to make the transi- 
tion diagram loop- free is unique. As an example, 
consider a TR diagram of Fig. [TIT a*). The elimi- 
nation of the edges e% : S2 — > S3 and 62:54—* 
S5 makes the diagram loop-free, as is illustrated 
in Fig. [TTT b). We can eliminate, say, the edges 
Sq — > S\ and S3 — > S4 to make the diagram loop- 
free. The elimination of two edges is necessary and 
sufficient to make the diagram loop- free in this case. 



uu 




Figure 11: An example with three loops, (a) origi- 
nal diagram, (b) modified diagram. 

Now, we recover the eliminated edges in the re- 
duced loop-free diagram one by one. At each re- 
covery step, at least one new loop is recovered and 
we select one such loop. These loops constitute ba- 
sis loops. In case of Fig. 1111 the recovery of edge 
ei : S2 — > S3 recovers the loop L x : So — > Si — > 
S2 ~^ S3 — > Sq- Then, recovery of e-i : S4 — > S5 
recovers L% So — ► S3 — > S4 — ► S5 — > Sq and 
L3 : So — > Si — > S2 — > S3 — > S4 — > S5 — > Sq. Thus, 
we can choose Li and L2 as basis loops. Li and L3 
can also be basis loops. In this respect, each edge 
selected to make the transition diagram loop-free 
corresponds to one basis loop. 

Now, let us consider an algebraic representation 
for edge removal. Consider an edge e : Sj — > Si. 
The probability flow of the edge e from Sj to Si is 
given by (|19|) . which is alternatively represented as 

Pij = v(i,j) T P, (51) 
where o~(i,j) is an (N + l)-vector given by 

( Qij, k=j + l 
Hi,j)} k = \ k = i + l (52) 

[ 0, otherwise. 



Note that pt appears at the [i + l)-th place instead 
of at the i-th place in the vector p because of the 
existence of po as its first element. Removal of e 
from the transition diagram amounts to removal of 
pij and pji = —pij from the [i + l)-th and the (j + 
l)-th components of Qp, respectively. Therefore, 
removal of e corresponds to removal of the matrix 

AQ:=-6(i,j)<j(i,j) T , (53) 

where S(i,j) is an (N + l)-dimensional vector given 

by 

( 1, k = i + 1 
[S(i,j)] k = l 1, k = j+l (54) 

The reduced transition matrix Q is given by 

Q' = Q - AQ. (55) 

Hence, from (|53|) and (|5~Tj) . we have 

Q'p = (Q + S(i,j)o-(i,j) T )p = PijS(i,j). 

Assume that the transition diagram contains I basis 
loops Li,Li2, ■ ■ ■ ,Li. We can assume that all these 
loops contain po without loss of generality. Now, 
choose a set of edges e k £ k = 1, 2, • • • ,1, such 
that the elimination of e\, e<2 ■ • ■ , e; makes the tran- 
sition diagram loop-free. We assume that the edge 
e/j connects the states S Mfc and S Vk in the direc- 
tion from S Mfc to S Uh , i.e., e k ■ S Mfc — » S Vk . Denote 
the probability flow associated with e k by p k , i.e., 
Pk = Pv kt i k - Then, from (J5jj), we get 

Pk = Qvk^kPfJ-k ~ QnkVkPvk 

= a(iy k ,p k ) T p. (56) 

The transition matrix Q' corresponding to the re- 
duced loop-free transition diagram is thus 

l 

Q' = Q + ^S(iy kl p k )a(u kl p k ) T . (57) 
fe=i 

Now we define a vector £(fk, Ph) satisfying 

Q'£,(vk,Pk) = 5{vk,Hk)> k = 1,2, • • • ,1. 
(K,^)o = 0. (58) 

We call £(j/fc,^tfc) the modifier of L k corresponding 
to the edge S^ k — > S 1 ^ , which is a generalization of 
£ introduced in the preceding section satisfying (|4"3")) 
or ([4*5]) . Since Q' corresponds to a loop- free transi- 
tion matirx, its stationary probability distribution 
p' is given by 

Pi = hoPo, i = 1,2, ■ ■ ■ ,N. (59) 
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where fio denotes the transition ratio from So to 
Si taken along a unique path connecting So to Si 
in the reduced loop-free transition diagram. From 
([57)1 and J55J), we get 

I 

Q'(i - ^2£,("k, Hk)v(vk, Vk) T ) = Q- 

k=l 

From Qp = 0, 

i 

fc=i 

From Q'p' = and the uniqueness of the stationary 
distribution, we have 

i 

fc=l 

Taking (|5ip into account, we can now write the sta- 
tionary probability distribution explicitly as 

i 

P = P + ^2pk£,(vk,Hk), (60) 

k=l 

where p' denotes the stationary probability distri- 
bution of the reduced loop-free transition diagram 
given by ([55)) . This generalizes (|44"]) . Equation 
([50)) clearly shows the importance of the probability 
flows reflecting the loop structure of the transition 
diagram. If the probability flows vanish for each 
loop, the stationary distribution is identical to that 
of reduced loop- free diagram given in ([317)) . The 
modifier £(i/k, p>k) represents the dependence of the 
stationary distribution on the probability flow. 

It remains to compute the probability flow 
pi,P2,--- , ph Premultiplication of ([6"0")) by 
a(v kl p k ) T yields 

( 

Ph ~ ^ cr(z/ fc , p k f £(V m , p. m )p m = v{vk,Pk) T p' '■ 

rn — 1 

(61) 

Now, from the definition of o~(i>k,pk) ([52p , we have 
a(v k ,Hk) T p' = t?^MfcP^ - Q^ k p'u„ 

= -Qn h v k r Vh0 (l - fo^rv^f^Po 

As is shown in Fig. 1121 the term ro Uk r Vk a k Tu k o corre- 
sponds to the product of the transition ratios along 
Lfc associated with the edge connecting S„ k and 
S Mfc , i.e., the Wegscheider product of L^. We de- 
note it by rjk- Equation (|61|) can now be rewritten 



as 

/ 

Pk+2_j ®kmPm = a*(l - 7?fe)po, fc = 1,2, • • • , I, 

m— 1 

(62) 

where 9km = -ff(^,Mi) T ?(vm,fi m ), a*: = 
—q^Vk^^Oi and 77^ is the Wegscheider product of 
-L/c . Equation (j62p implies that the probability flows 
can be computed by solving only a linear equation 
of I unknowns, instead of (N + 1) unknowns of SSE 
© or (gj). We call ||52J) the /Zow equation. 




Figure 12: Loop Elimination by Removal of an 
Edge. 

Now, let us summarize the procedure of comput- 
ing the stationary probability distribution. 

Step 1. Find the minimum number of edges e k '■ 
S/j, k — » S Vk , k — 1, 2, • • • , I such that their elimina- 
tion makes the transition diagram loop-free. 

Step 2. Compute the stationary distribution p' via 
([59)) for the reduced loop-free diagram. 

Step 3. For each eliminated edge e^, associate a 
loop Lk, which is recovered with the edge. 

Step 4. For e*, compute the modifier £(i>k, Mfc) that 
satisfies ([55)) . 

Step 5. Solve the flow equation ([6"2")) to obtain the 
probability flows pk, k = 1, 2, • • • , I. 

Step 6. Compute p from (|60p . 

Step 7. Compute the exact probability distribu- 
tion taking the normalization constraint ([5]) into 
account. 

The computation of the modifier (Step [4} is dis- 
cussed in Appendix. 

Example 7. Consider the transition diagram 
shown in Fig. [TITa) with 6 states and 3 loops. To 
eliminate the loops, two edges, e\ : S2 — > S3 and 
62 : S4 — > S5, are removed. The modified transition 
diagram is shown in Fig. Illf b). which is loop- free. 
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We choose L\ : So — > Si — > S2 — > S3 — > So, which 
is created by recovering ei in the modified diagram 
and L2 : Sq — > S3 — > 54 — > S5 — > So which is cre- 
ated by recovering e2 in the modified diagram as 
basis loops. In terms of the notations introduced 
in this section, [i\ = 2, v\ = 3, p,2 = 4, v<i = 5, 
Pi = P32 and P2 = /O54. The transition matrices Q 
and Q' for the original diagram and the modified 
diagram are given by 



rji and 772 are Wegscheider products corresponding 
to Li and L2, respectively, and are given respec- 
tively: 



Vi = ^10^21^32^03, V2 = r 30 r4 3 r M r 05 . 



Q = 



Q' = 



D a 


9oi 





903 





905 


9io 


Di 


912 














921 


D 2 


1923 i 








930 





1932; 


D 3 


934 














943 


£>4 


|945| 


950 











|954| 


D 5 _ 


".Do 


901 





903 





905 




9io 


Di 


912 
















921 


D' 2 













930 








D's 


934 
















943 


D\ 







.950 














D' 5 _ 





-1- 930 +950), 

D 3 = -(903 

-(905 + 945 ] 



f 923 

D' = 



-(901 - 
- 943), 
-912, 



-921] 

D 4 -- 
D' = 



where D — — (qio 
D 2 = -(912 + 932), 
-(934 + 954), A; = 
-(903 + 943), D' 4 = -q 3i , D' b = -q 05 . The elements 
of Q encircled by dashed boxes are eliminated in Q' . 
The modifiers, £1 corresponding to Li for ei and £2 
corresponding to L2 for e2, are given as 





I/901 
-^21/901 - I/912 

1/903 

?W903 









-1/903 
-^•43/903 - 1/934 
1/905 

(63) 



The detail of the above calculation is given in Ap- 
pendix. 

According to (|62|) . the flow equation is given by 
(1 + 6 n )pi + 812P2 = «i(l - m) 

82lPl + (1 + #22)P2 = "2(1 - V2) 

where 9^ = —of^j, i,j — 1,2, with 



a(3,2f 
a(5,4f 



[0 g 3 2 -923 0] 
[0 q 54 -,745 0] 



and 



Example 8. (Enzyme with three effectors) 
Consider a more complex block transition diagram 
shown in Fig. I13f a). This transition diagram de- 
scribes an allosteric enzyme with three different ef- 
fectors or an operon with three transcription fac- 
tors. A description like (|16p for this example is very 
complicated and hence, is omitted. Now, we choose 
four edges e\ : Si — > S 2 , e 2 : S 4 — > S 5 , e 3 : S 6 — » S 7 
and e4 : S7 — * Sg to make the diagram loop-free. 
The reduced diagram is shown in Fig. [T37 b). A basis 
set of loops is composed of the following four loops, 
Li ■ So — > Si — > 52 — * S3 — > So, L2 : So — * S3 — » 
S4 — » S5 — > So, £3 : So — > Si — * 5*6 — » 5*7 — > So 
and L4 : So — > S7 — > Ss — * S5 — > So, as is shown 
in Fig. I13f a). The stationary probability distribu- 
tion corresponding to the reduced loop- free diagram 
Fig. rfBT b) is easily calculated according to (|59|) . i.e., 



Pi 

Pa 
P' e 



rioPo, P 2 = r 2oPo = r2 3 r 30 po, p 3 = r 3a po, 
fioPo = r i3 r 3a po, p' b = r 50 po, 
feoPo = reirioPo, r' 7 = r 70 p , 
rsoPo = r S5 r 50 p . 



The a- vectors defined by (|52[) in this case are o\ = 
a{2, 1), <7 2 = <r(5, 4), a 3 = a{7, 6), <t 4 = a{8, 7), i.e., 



<Tl 





















921 

















-912 








































954 


, 0-3 = 





, cr 4 = 










-945 






















976 

















-967 




987 

















-978 



o-i 



"923?"30, OL2 



-945?"50- 



The modifiers £j, i = 1, 2, 3, 4 for L\, L2, L 3 , L 4 are 
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computed in Appendix; they are 



6 = 



£3 = 











-1/901 







?W<703 + V932 




^3/903 


1/903 




- 1/903 


^43/903 


, 6 = 




-^43/903 - 1/934 







1/905 






















?W905 










- 1/901 


























, & = 












I/905 


-r el /q 01 - l/q w 







1/907 




- 1/907 







735/905 + 1/954 



(64) 



The flow equation 



l + 6»] 

0-21 
031 



Qfl(l - 

«2(i - m) 
c*3(i - m) 
«4(i - m) 



is given by 



U\2 


Wl3 


l+9 22 








1 + 6*33 


042 


^43 





024 
#34 
+■ #44 





"pl" 




P2 




P3 




P4. 



(65) 



where 



= 1,2,3,4, j = 1,2,3,4, 

«1 = -912^20, a 2 = -945^50, a 3 = -967^70, 

04 = —978^80, and 77.;, i = 1,2, 3, 4, are Wegscheider 
products for Li, and are given by r/i — r^r^rzirio, 

m = ^05^54^43^30, % = r Q7 r 76 r el r w , and Tf 4 = 

r 05 r 58 r 87'''70- The anti-diagonal elements 9u, #23, 
632 and 04i vanish from the loop structure of Fig. [13] 



7 Specification of transition 
probabilities 

In order to obtain a useful representation of the ac- 
tivity of the biological regulator discussed so far, we 
must give a way to connect transition probabilities 
qij with more workable parameters with physical 
and/or biological meanings. The state transition 
caused by binding a BF to a BS and that caused 
by releasing a BF from a BS are essentially differ- 
ent in their nature. The former transition is usually 
proportional to the concentration of the BF to be 
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Figure 13: An Example of Transition Diagram, (a) 
Original Diagram, (b) Reduced Loop-free Diagram. 
Transition probabilities are omitted. 



bound, while the latter does not depend on concen- 
tration of BFs but possibly depends on the occupa- 
tion pattern of other BSs. Therefore, it is natural to 
assume the following characterization of transition 
probabilities: 



qij 



XijUk, if Sj — » Si is the binding of Uk, 



for each j and i € A.j , where a^- 



if it is a releasing of a BF 

(66) 

denotes a coeffi- 
cient and Uk the concentration of Uk ■ The transition 
probabilities used in Example 1 (equation (|10p ) and 
these in Example 2 (equation HH)) are special cases 

of HMD- 

From the definition of TR in (12"4"1) , we have 



_ / PijUk, 



A 



if Sj — ► Si is the binding of Uk, 
if it is the releasing of Uk- 

(67) 

where ftij — ctij/aji. We can derive a concrete 
form for the overall activity of our biological regu- 
lator based on (|66|) , for the case that the Wegschei- 
der condition holds, i.e., the stationary distribution 
given by (|3Uj) . 

Let Si be a state where I sites are occupied by U^ , 
Ui 2 , ■ ■ • , Ui t and the remaining n — I sites are empty. 
Then, according to (|28| and (|66| . the probability 



Pi IS = CliU^Uiz 



where a,- denotes the 
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product of flijS associated with the binding of Z7 ii; 
Ui 2 , ■ ■ • , Ui t . Note that the binding and unbinding 
of other BFs cannot take place because the path 
connecting So to Si does not contain any loop. Now, 
the normalization constraint ([5]) yields the overall 
activity 7 defined in ([6|) as 

_ 7o + £iIi7»QWi^2 ■■■ u h (m) 
1 -\- ij i= ^aiUi 1 Ui 2 ■ ■ ■ Ui l 

This is a familiar form for describing operon reg- 
ulation, which appears frequently in the literature 
e.g., [30] [34] 00]. For the Lac operon of Fig. 
the transition probabilities given by (|66[) implies 
9io = a w ui, q i = a i, q 2 \ = a 2X u 2 , q\ 2 = (X12, 
<?30 = a 30 u 2 , qo3 = «03, <?40 = "40^3 and q 04 = a aA . 
Since only the states bound by U\ can initiate tran- 
scription, 70 = 73 = 74 = in ([68]) . Thus, we have 

_ liKiui + j 2 KiK 2 um 2 

7 ~~ 1 + K lUl + K x K 2 uiu 2 + K 3 u 2 + K 4 u 3 ' 

where K\ = a w /a i, K 2 = a. 2 \ja 12l K 3 = a 30 /a 03 
and K4 = Q!4o/ao4- 

The general form (|4Tj) , as well as its specification 
([68]) . enables us to derive a variety of complicated 
formulae representing enzymic actions (e.g., [40] ) 
almost immediately. Equation (|68[) also generalizes 
the classical MWC model [35] to heterotropic cases. 

It is important to notice that the Wegscheider 
product is constant and does not depend on concen- 
tration of any transcription factor given the spec- 
ification ([66]) . To see this, notice that if the loop 
contains an edge associated with the binding of Uk, 
it must contain an edge associated with unbinding 
of Uk, because the loop must recover the starting 
state along the path. If the binding transition con- 
tains a TR with Uk, then it must contain TR as- 
sociated with unbinding of Uk which contains , 
so that they cancel out in the Wegscheider product. 
Thus, we have shown that the Wegscheider product 
is always constant and does not contain the concen- 
trations of BF. 

8 Conclusions 

The metabolic process is controlled by enzymes 
whose expressions are controlled by genetic regu- 
lations. On the other hand, genetic regulation is 
controlled by the products of metabolism that de- 
termine the cell state. In this sense, the genetic and 
metabolic regulations are closely linked together to 
form a huge and complex network of intracellu- 
lar regulations. It has been desirable to establish 
a common framework for quantitatively describing 



these regulations in a unified way. For that purpose, 
we used the analogy between genetic and metabolic 
regulations shown in Table [T] The core of the anal- 
ogy is that the actual computations of the con- 
trol action are performed through molecular inter- 
actions at the regulatory sites between the sites to 
be bound and the factors to bind or to dissociate. 

We formulated a finite Markov process describing 
both the genetic and metabolic regulations. The 
most salient feature of this Markov process is its 
reciprocity in transition probability; that is, for 
each pair of state, if the transition probability from 
one state to another is positive, the opposite tran- 
sition probability is also positive. This property re- 
flects the reversibility of the molecular interactions 
we are dealing with, and it is the source of many in- 
teresting properties of the Markov process we have 
formulated. The parameters that determine the ac- 
tion of regulations is represented as the average rate 
of the stationary probability distribution of that 
Markov process, based on the assumption that the 
dynamics of molecular interactions that compute 
the control actions are much faster than the reac- 
tion they are regulating. This assumption is analo- 
gous to the fast equilibrium assumption used in de- 
riving the classical Michaelis-Menten equation [40] . 
Actually, our approach can derive a stochastic ver- 
sion of Michaelis-Menten formula which suggests its 
legitimacy. Moreover, we can derive a quantitative 
estimate of how precise the Michaelis-Menten equa- 
tion is by supplying the variance of the stationary 
distribution. 

We introduced a new notion of probability flow 
associated with each edge of the transition diagram 
to represent the loop structure of the transition di- 
agram. The state transition diagram is nothing but 
a representation of the conservation of probability 
flows at each node. Based on this observation, we 
derived many interesting properties of the station- 
ary probability distributions. 

We proved a simple result that the probability 
flows vanish if the graph has no loop. This imme- 
diately implies that the detailed balance holds for 
pairs of adjacent states. This was generalized to in- 
clude pairs of states which are not necessarily adja- 
cent to each other. A very simple graphical method 
of computing the stationary distribution was de- 
rived. 

We derived a condition which guarantees the de- 
tailed balance even if the transition diagram has 
a loop. This condition is not new; Wegscheider, 
a German chemist, discovered this condition a cen- 
tury ago. The condition, which we call the Wegshei- 
der condition in this paper, requires that the prod- 
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uct of all the transition ratios around a loop be 
unity. We combined this condition with the proba- 
bilistic flow and directly showed that the probability 
flows vanish under it. This result again gives a prob- 
abilistic interpretation of this classical result. The 
detailed balance has been accepted as an obvious 
fact in literature of chemical kinetics and enzymol- 
ogy, For instance, a famous book of Segel assumed 
this condition without even mentioning Wegschei- 
der's name [40] . We gave a stochastic interpretation 
of this condition. 

We derived a simple method of computing the 
stationary probability distribution for general cases 
with multiple loops. We showed that if the Wegshei- 
der condition holds for a loop, we can eliminate any 
of the edges contained in the loop without chang- 
ing the stationary probability distribution. In other 
words, if the Wegscheider condition holds for a loop, 
we can consider a transition diagram with that loop 
by removing an arbitrary edge within that loop. 
The stationary probability distribution of the re- 
duced diagram is identical to that of the original 
diagram. 

We derived a purely graphical method of comput- 
ing the stationary probability distribution based on 
the notions of probability flow and modifiers. We 
derived a simple equation for computing the proba- 
bility flow. Our method dramatically simplifies the 
classical King and Altman method [24] . 

Our result suggests that there are two kinds of 
stationary probability distributions: the one which 
satisfies the detailed balance, and one which does 
not. The former class is characterized by zero prob- 
ability flow, the latter by non-zero probability flows. 
In the latter case, the probability distribution is 
fixed, but continuous flows of probability exist in 
loops, which perhaps reflects a sort of thermal irre- 
versibility of molecular interactions. Exploitation of 
the bio-chemical characterization of the probability 
flow is an interesting issue of theoretical biology. 

The theoretical framework developed in this pa- 
per is based on the observation that the intracellular 
regulations are embedded in a homogeneous com- 
putational medium of molecular interactions. This 
view is not entirely new, but no serious attempt has 
been made so far to mathematically formulate it, 
within the best of our knowledge. A finite Markov 
process model proposed in this paper captures some 
essential features of the computational medium and 
explains the versatility, evolvability and flexibility 
of the intracellular regulations and various signal 
transductions. It offers a potential capability to 
deal with systems in which metabolism and gene 
expressions are linked and integrated together. We 



are now exploiting an analytical tool for investigat- 
ing such systems with greater complexity based on 
our framework presented in this paper. 

References 

[1] Ackers, O.K., Johnson, A.D., Shea, M.A., 
1982. Quantitative model for gene regulation 
by A phage repressor. Proc. Nat. Acad. Sci. 79, 
1129-1133. 

[2] Beckett, D., 2005. Multilevel regulation of 
protein-protein interactions in biological cir- 
cuitry. Phys. Biol. 2, S67-73. 

[3] Berg, O.G., von Hippcl, P.H., 1987. Selection 
of DNA binding sites by regulatory proteins. 
Statistical-mechanical theory and application 
to operators and promoters. J. Mol. Biol. 193, 
723-750. 

[4] Bintu, L., Buchler, N.E., Garcia, H.G., Ger- 
land, U., Hwa, T., Kondev, J., Phillips, R., 
2005. Transcriptional regulation by the num- 
bers: models. Curr. Opin. Genet. Dev. 15, 116- 
24. 

[5] Bliss, R.D., Painter, P.R., Marr, A.G., 1982. 
Role of feedback inhibition in stabilizing the 
classical operon. J. Theor. Biol. 97, 177-193. 

[6] Bolouri, H., Davidson, E.H., 2002. Modeling 
DNA sequence-based cis-regulatory gene net- 
works. Dev. Biol. 246, 2-13. 

[7] Chen, K.C., Csikasz-Nagy, A., Gyorffy, B., Val, 
J., Novak, B., Tyson, J.J., 2000. Kinetic anal- 
ysis of a molecular model of the budding yeast 
cell cycle. Mol. Biol. Cell 11, 369-91. 

[8] Cherry, J.L., Adler, F.R., 2000. How to make a 
biological switch. J. Theor. Biol. 203, 117-133. 

[9] Coppey, M., Benichou, O., Voituriez, R., 
Moreau, M., 2004. Kinetics of target site lo- 
calization of a protein on DNA: a stochastic 
approach. Biophys. J. 87, 1640-1649. 

[10] Dekel, E., Alon, U., 2005. Optimality and evo- 
lutionary tuning of the expression level of a 
protein. Nature 436, 588-592. 

[11] Endy, D., Brent, R., 2001. Modelling cellular 
behaviour. Nature 409, 391-395. 

[12] Fell, D., 1997. Understanding the Control of 
Metabolism. Portland Press, London. 



17 



[13] Freire E., 2000. Can allosteric regulation bo 
predicted from structure? Proc. Natl. Acad. 
Sci. USA. 97, 11680-11682. 

[14] Gardner, T.S., Cantor, C.R., Collins, J.J., 
2000. Construction of a genetic toggle switch 
in Escherichia coli. Nature 403, 339-342. 

[15] Gibson, K.M., Gircn, J.A., Head, M.S., 1997. 
A new class of models for computing receptor- 
ligand binding affinities. Chem. Biol. 4, 87-92. 

[16] Gillespie, D.T., 1977. Stochastic simulation of 
coupled chemical reactions. J. Phys. Chem. 81, 
2341-2361. 

[17] Griffith, J.S., 1968. Mathematics of cellular 
control processes. I. Negative feedback to one 
gene, II. Positive feedback to one gene. J. 
Theor. Biol. 20, 202-208, 209-16. 

[18] Gunasekaran, K., Ma, B., Nussinov, R., 2004. 
Is allostery an intrinsic property of all dynamic 
proteins? Proteins 57, 433-43. 

[19] Heinrich, H., Schuster, S., 1996. The regulation 
of Cellular Systems. Chapman and Hall, New 
York. 

[20] Istrail, S., Davidson, E.H., 2005. Logic func- 
tions of the genomic cis-regulatory code. Proc. 
Natl. Acad. Sci. USA. 102, 4954-4959. 

[21] Kaern, M., Elston, T.C., Blake, W.J., Collins, 
J.J., 2005. Stochasticity in gene expression: 
from theories to phenotypes. Nat. Rev. Genet- 
ics 6, 451-464. 

[22] Karmarkar, R., Bose, I., 2004. Graded and bi- 
nary responses in stochastic gene expressions. 
Phys. Biol. 1, 197-204. 

[23] Killer, A.D., 1995. Model genetic circuit en- 
coding autoregulatory transcription factors. J. 
Theor. Biol. 172, 169-185. 

[24] King, E.L., Altman, C, 1956. A schematic 
method of deriving the rate laws for enzyme- 
catalyzed reactions. J. Phys. Chem. 60, 1375- 
1378. 

[25] Kitano, H., 2004. Biological robustness. Nat. 
Rev. Genet. 5, 826-837. 

[26] Lloyd, G., Landini, P., Busby, S., 2001. Acti- 
vation and repression of transcription initiation 
in bacteria. Essays Biochem. 37, 17-31. 



[27] Mackey, M.C., Santillan, M., Yildirim, N., 
2004. Modeling opcron dynamics: the trypto- 
phan and lactose operons as paradigms. C.R. 
Biol. 327, 211-24. 

[28] McAdams, H.H., Arkin, A., 1997. Stochas- 
tic mechanisms in gene expression. Proc. Nat. 
Acad. Sci. 94, 814-819. 

[29] Ming, D., Wall, M.E., 2005. Quantifying al- 
losteric effects in proteins. Proteins 59, 697- 
707. 

[30] Mochizuki, A., 2005. An analytical study of 
the number of steady states in gene regulatory 
networks. J. Theor. Biol. 236, 291-310. 

[31] Monod, J., 1972. Chance and Necessity, Paper- 
back. 

[32] Monod, J., Wyman, J., Changeux, J.P., 1965. 
On the nature of allosteric transitions: a plau- 
sible model. J. Mol. Biol. 12, 88-118. 

[33] Munsky, B., Khammash, M., 2006. The finite 
state projection algorithm for the solution of 
the chemical master equation. J. Chem. Phys. 
124, 044104 

[34] Ozbudak, E.M., Thattal, M., Lim, H.M., 
Shralman, B.I., von Oudenaarden, A., 2004. 
Multistability in the lactose utilization network 
of Escherichia coli. Nature 427, 737-740 

[35] Pan, H., Lee, J.C., Hilser, V.J., 2000. Bind- 
ing sites in Escherichia coli dihydrofolate re- 
ductase communicate by modulating the con- 
formational ensemble. Proc. Natl. Acad. Sci. 
USA 97, 12020-12025. 

[36] Perutz, M.F., 1989. Mechanisms of coopera- 
tivity and allosteric regulation in proteins. Q. 
Rev. Biophys. 22, 139-237. 

[37] Ptashne, M., 1992. A Genetic Switch, Cell & 
Blackwell Scientific, Cambridge, MA. 

[38] Santillan, M., Mackey, M., 2004. Influence of 
catabolite repression and inducer exclusion on 
the bistable behavior of the lac operon. Bio- 
phys. J. 86, 75-84. 

[39] Savageau, M.A., 1985. A theory of alterna- 
tive designs for biochemical control systems. 
Biomcd. Biochim. Acta. 44, 875-880. 

[40] Segel, I.H., 1993. Enzyme Kinetics: Behaviour 
and analysis of rapid equilibrium and steady- 
state enzyme systems, John Wiley and Sons, 
INC., New York. 



18 



[41] Smolen, P., Boxter, D.A., Byrne, J.H., 
2000. Modeling transcriptional control in gene 
networks-methods, recent results, and future 
directions. Bull. Math. Biol. 62, 247-292. 

[42] Tanaka, R.J., Okano, H., Kimura, H., 2006. 
Mathematical description of gene regulatory 
units. Biophys. J. 91, 1235-1247. 

[43] Thattai M., van Oudenaarden, A., 2001. In- 
trinsic noise in gene regulatory networks. Proc. 
Natl. Acad. Sci. USA 98, 8614-8619. 

[44] Thomas, R., Thieffry, D. Kauffman, M., 1995. 
Dynamical behaviour of biological regulatory 
network - I. Biological role of feedback loops 
and practical use of the concept of the loop- 
characteristics state. Bull. Math. Biol. 57, 247- 
276. 

[45] Tyson, J.J., 1978. The dynamics of feedback 
control circuits in biochemical pathways. Prog. 
Theol. Biol. 5, 1-61. 

[46] van Kampen, N.G., 1992. Stochastic Process in 
Physics and Chemistry. Elsevier Science, Am- 
sterdam. 

[47] Wolf, D.M., Eeckman, F.H., 1998. On the re- 
lationship between genomic regulatory element 
organization and gene regulatory dynamics. J. 
Theor. Biol. 195, 167-186. 

[48] Wong, P., Gladney, S., Keasling, J.D., 1997. 
Mathematical model of the lac operon: in- 
ducer exclusion, catabolite repression, and di- 
auxic growth on glucose and lactose. Biotech- 
nol. Prog. 13, 132-143. 

[49] Wyman, J., 1968. Linked functions and recip- 
rocal effects in hemoglobin: A second look. 
Adv. Protein Chem. 19, 223-286. 

[50] Yanofsky, C, 1992. Transcriptional regulation: 
Elegance in design and discovery. Transcription 
Regulation I. Cold Spring Harbor Lab. Press, 
pp.3-24. 

[51] Yi, T.M., Huang, Y., Simon, M.L, Doyle, J., 
2000. Robust perfect adaptation in bacterial 
chemotaxis through integral feedback control. 
Proc. Natl. Acad. Sci. USA. 97, 4649-4653. 

[52] Zaslaver, A., Mayo, A.E., Rosenberg, R., 
Bashkin, P., Sberro, H., Tsalyuk, M., Surette, 
M.G., Alon, U., 2004 Just-in-time transcrip- 
tion program in metabolic pathways. Nat. 
Genet.36, 486-491. 



Appendix 

Computation of the modifier can be done through 
the simple recursions presented in Section [4] The 
modifier corresponding to the loop Lk is an (N+l)- 
dimensional vector — £,(vk,^k) satisfying ([58]) . 
For the case that the transition diagram itself is 
a loop, we have seen that the recursions formulae 
(|46|) and ([47]) give a vector £ that satisfies (|48|). 
which corresponds to the special case [i^ — v, vu — 
v+1. The computation of the modifier for a general 
loop is essentially reduced to (|46f (|47|) . To avoid 
notational complication, we assume that Lk is the 
loop So — > Si — > 52 — > • • • — > — > So- Since 
efc : S^ lk — > S Uk is in the loop, we assume fik = v, 
Vk = v + 1, v < 7T — 1. Then, the components 
of the modifier corresponding to the states inside 
the loop Lk are given by (l46|) and (j47|) . Thus, we 
have a procedure for computing the components of 
the modifier corresponding to the states inside the 
loop. The components of outside the loop Lk are 
computed as follows: let Sj be a state outside Lk- 
Since the reduced diagram is loop-free, there exists 
a unique path connecting Sj to So- If the path does 
not have a common state with Lk except So, let the 
component of corresponding to Sj be zero. If the 
path has some states in common with Lk, there is 
a state Sfc, nearest Sj in Lk- Then, the component 
of corresponding to Sj is given by 

(&)i=* ! Jfc 1 & J ( Al ) 

where fjkj is the TR from the state Sfc, to Sj and 
. is the component of the modifier corresponding 
to Sk which has already been computed. The jus- 
tification of (Al) is based on the remark at the end 
of Section [3] (equation (|3~Tj) ). 

Example Al. We compute £i and £2 in Fig. [TT1 
For L±, the forward recursion corresponding to (|46[) 
is given by 

60 = 



£l2 = r2l£ll = -T21 

Q12 9oi 912 

and the backward recursion corresponding to (|47[) 
is given by 



From the graph of Fig. [TlTb') , the path connecting 
S4 to So meets L\ at S3, while the path connecting 
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Sg to Sq directly reaches So- Hence, we have 
64=^43^13 = —, £i5 = 0. 

903 

For L2, the forward recursion is given by 
60 =0 
£23 = 

903 

£24 = ^43 £23 



7-43 1_ 

903 934 



and the backward recursion is given by 
£25 = — • 

905 

For the states outside L 2 , S2 and Si are directly 
connected to So without meeting L 2 . Hence, £22 — 
£21 = 0. Thus, ([55]) has been confirmed. 

Example A2. We compute £1, £2, £3 and £4 corre- 
sponding to the loops L±, L 2 , Li and L4 of Example 
8. 

For L±, the forward recursion (Sq — > Si) gives 

£11 = -I/901, 
and the backward recursion (So ■ — > S3 — > S2) gives 

63 = 1/903, 

£12 = r 23 £,13 + 1/932- 

The probabilities of the states outside L\ are given 

by 

£14 = r 43 £i 3 , 

iu = f i0 £w = 0, i = 5, 6, 7, 8. 

For L 2 , the forward recursion (Sq — > S3 — > S4) 
gives 

£23 = ~ I/903, 
£24 = r 43 £ 2 3 - I/934, 

and the backward recursion (So —> S5) gives 

65 = 1/905- 

The probabilities outside L 2 are given by 

£22 = ^23£23, 
^28 = r 85£25, 

6i = ho&a = 0, i = 1,6,7. 



For L3, the forward recursion (5*0 — > Si — ► Sg) 
gives 

£31 = -I/901, 

62 = rn £31 - I/914, 

and the backward recursion (So — ► S7) gives 

£37 = I/907 
. The probabilities outside L3 is given by 
£3* = fio6o = 0, » = 2, 3, 4, 5, 8 

For L4, the forward recursion (So — ^ S7) gives 

£47 = - 1/907 
, and the backward recursion (So — * S5 — ► Sg) gives 

£45 = I/905, 

£48 = r S5 £45 + 1/954 • 

The probabilities outside L4 is given by 

£41 = no&o = 0, i = 1, 2, 3,4, 6. 

The above scheme verifies (IMl) . 
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